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In  recent  years,  vulnerability  against  high-velocity  impact  loads  has  become  an  increasingly  critical  issue 
in  the  design  of  composite  aerospace  structures.  The  effects  of  Hydrodynamic  Ram  (HRAM),  a  phenom¬ 
enon  that  occurs  when  a  high-energy  object  penetrates  a  fluid-filled  container,  are  of  particular  concern 
in  the  design  of  wing  fuel  tanks  for  aircraft  because  it  has  been  identified  as  one  of  the  important  factors 
in  aircraft  vulnerability.  The  projectile  transfers  its  momentum  and  kinetic  energy  through  the  fluid  to 
the  surrounding  structure,  increasing  the  risk  of  catastrophic  failure.  In  the  present  paper,  the  commercial 
finite-element  code  ABAQUS/Explicit  has  been  used  to  simulate  an  HRAM  event  due  to  the  impact  of  a 
steel  spherical  projectile  into  a  water-filled  woven  CFRP  square  tube.  In  order  to  simulate  the  fluid-struc¬ 
ture  interaction,  the  Coupled  Eulerian  Lagrangian  (CEL)  approach  is  used.  Experimental  tests  which  indi¬ 
cate  the  pressure  at  different  points  of  the  fluid,  strains  of  the  walls  and  cavity  evolution  for  different 
impact  velocities  are  compared  with  the  numerical  results  in  order  to  assess  the  validity  and  accuracy 
of  CEL  technique  in  reproducing  such  a  complex  phenomenon.  Also,  several  numerical  impacts  at  differ¬ 
ent  initial  projectile  velocities  are  performed  to  study  its  influence  in  the  HRAM  phenomenon. 

©  2012  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Nowadays  the  aeronautical  and  aerospace  industries  are  contin¬ 
uously  increasing  the  usage  of  laminated  composite  structures  to 
diminish  the  fuel  consumption.  The  choice  of  this  kind  of  materials 
is  due  to  the  high  strength-to-weight  and  stiffness-to-weight  ratios 
as  well  as  their  anisotropic  behavior.  Those  special  characteristics 
allow  to  optimize  designs  and  fulfil  the  strict  requirements  of  the 
mentioned  industries  with  a  reduction  of  the  total  mass  of  the 
structure  and  hence  saving  fuel.  The  most  used  composite  materi¬ 
als  for  structural  applications  in  these  sectors  are  carbon  fiber  rein¬ 
forced  plastics  (CFRPs),  commonly  manufactured  with  an  epoxy 
matrix  that  combines  good  mechanical  properties,  high  resistance 
to  corrosion  and  fatigue,  and  low  density  (p  =  1500  kg/m3).  The 
reduction  of  raw  material  costs,  the  development  of  automation 
of  manufacturing  processes,  and  the  growing  experience  in  design 
technology  have  increased  the  CFRP  applications  in  commercial 
aircraft  [1],  such  as  the  fuselage  and  wings,  so  that  in  the  last  de¬ 
signs  of  aeronautical  structures  this  kind  of  materials  constitutes 
more  than  the  50%  (in  terms  of  weight).  The  successful  usage  of 
these  materials  in  primary  structures  depends  on  understanding 
their  response  to  a  wide  range  of  impact  loadings.  The  study  of 
the  behavior  of  composite  laminates  under  high-velocity  impact 
has  received  attention  since  the  1970s,  when  only  military 
research  agencies  where  concerned  about  this  subject.  In  the 
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1980s,  investigators  from  civil  organizations  began  to  publish 
articles  in  which  the  breakage  mechanism  of  such  materials  under 
ballistic  impact  was  investigated  primarily  from  an  experimental 
approach  [2-4].  Later,  both  analytical  [5-11]  and  numerical  [12- 
17]  approaches  were  used  to  predict  the  energy  absorbed  by  the 
laminate  in  a  penetrating  impact  and  the  damaged  area. 

Vulnerability  studies  of  CFRP  aerospace  structures  are  becom¬ 
ing  an  issue  of  great  importance  in  the  design  of  any  aircraft  [8]. 
These  structures  may  suffer  high  velocity  impact  loads  due  to  bird 
strikes  [18]  or  hailstones  [19],  especially  dangerous  because  of 
their  high  possibility  of  occurrence  and  their  disastrous  conse¬ 
quences.  Runway  debris  may  impact  the  underside  of  the  wing 
structures  [20]  causing  Hydrodynamic  Ram  effects  in  the  fuel 
tanks,  which  are  considered  one  of  the  most  important  factors  in 
aircraft  vulnerability. 

The  Hydrodynamic  Ram  (HRAM)  phenomenon  appears  when  an 
object  with  high  kinetic  energy  penetrates  a  fluid-filled  tank  and 
transfers  its  kinetic  energy  through  the  fluid  to  the  surrounding 
structure,  increasing  the  risk  of  catastrophic  failure  and  excessive 
structural  damage.  HRAM  is  particularly  dangerous  for  aircrafts 
with  lightweight  designs,  because  the  structural  resistance  of  their 
integral  fuel  tanks  cannot  be  improved  by  strengthening  the  air¬ 
frame;  strengthening  the  frame  would  counteract  the  require¬ 
ments  of  a  lightweight  design.  Vulnerability  to  HRAM  has  been 
usually  related  to  military  aircraft,  but  commercial  airplanes  are 
not  exempt  of  its  effect.  In  1990  the  Federal  Aviation  Administra¬ 
tion  (FAA)  established  the  Aircraft  Catastrophic  Failure  Prevention 
Research  Program  in  which  the  analysis  of  the  effects  of  an 
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(a)  FiberFailure  (b)  Matrix  in  plane  failure 

Fig.  2.  Stress  components  that  participates  in  each 


(c)  Matrix  Crushing  and  de¬ 
lamination 

failure  mode. 


uncontained  turbine  engine  fragment  penetrating  commercial  air¬ 
craft  fuel  tanks  [21  ]  was  carried  on.  An  example  of  the  importance 
of  the  HRAM  phenomenon  is  the  Concorde  accident  that  occurred 
in  2000.  The  final  investigation  report  revealed  that  the  HRAM  had 
played  a  significant  role  in  the  aircraft  failure. 

The  study  of  the  HRAM  phenomenon  is  not  only  important  for 
the  aircraft  industry.  High  velocity  impacts  on  fluid  filled  contain¬ 
ers  are  of  great  interest  for  different  industrial  fields  such  as  safety 
of  industrial  facilities  or  road  haulage,  where  containers  are  com¬ 
monly  used  to  store  fuel  or  dangerous  products.  In  those  cases, 
an  impact  in  the  vessel  may  produce  the  failure  of  the  tank  and 


could  result  in  a  serious  consequences  on  the  environment  or  even 
toxic  and  flammability  effects  [22-24]. 

Hydrodynamic  Ram  consists  of  four  principal  stages:  shock, 
drag,  cavitation  and  exit.  Each  stage  contributes  to  structural  dam¬ 
age  through  a  different  mechanism  and  to  a  different  extent.  When 
the  projectile  penetrates  the  wall  of  the  fluid  filled  structure,  the 
impact  energy  is  transferred  to  the  fluid  and  generates  a  high-pres¬ 
sure  hemispherical  shock  wave.  This  leads  to  damage  primarily  in 
the  vicinity  of  the  impact  position.  During  the  drag  phase,  the  pro¬ 
jectile  travels  through  the  fluid,  and  its  kinetic  energy  is  partially 
transformed  into  fluid  motion  as  the  projectile  is  slowed  by  viscous 
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Table  1 

Properties  of  woven  carbon/epoxy  laminate. 


Elastic  properties 

e,=e2 

e3 

V12 

68  GPa 

10  GPa 

0.22 

Vis  =  v23 

Ci  2 

G23  =  G13 

0.49 

5  GPa 

4.5  GPa 

Strength  properties 

X 

ii 

ii 

nX 

ii 

Zc 

Zr 

880  MPa 

340  MPa 

96  MPa 

$12 

$13 

S 23 

84  MPa 

120  MPa 

120  MPa 

Critical  strain 

£l  =  £2 

S3 

£12  =  £23  =  £13 

0.025 

0.05 

0.1 

drag.  The  displacement  of  the  fluid  from  the  projectile  path  gener¬ 
ates  a  radial  pressure  field.  In  contrast  to  the  pressure  field  devel¬ 
oped  during  the  shock  phase,  the  fluid  is  accelerated  gradually 
rather  than  impulsively.  This  causes  less  intense  peak  pressures, 
but  they  are  of  greater  temporal  extent.  The  displacement  of  fluid 
during  the  drag  stage  forms  a  cavity  behind  the  projectile.  The  sub¬ 
sequent  expansion  and  collapse  (oscillations)  of  the  cavity  is 
known  as  the  cavitation  stage.  The  oscillations  of  the  cavity  can 
cause  significant  pressure  pulses.  The  final  stage  of  Hydrodynamic 
Ram  occurs  when  the  projectile  exits  the  container.  In  contrast  to 
the  perforation  of  the  front  wall,  the  exit  of  the  projectile  occurs 
through  a  pre-stressed  wall.  The  pre-stress  is  caused  by  the  initial 
shock  stage  and  the  subsequent  loading  by  the  fluid. 

Simulation  of  HRAM  events  has  been  attempted,  with  more  or 
less  success,  for  over  30  years.  The  first  methods  were  based  on 
the  use  of  the  Piston  Theory  [25,26]  and  the  Variable  Image  Method 
[27]  for  the  fluid-structure  interaction.  Once  these  methods  were 
proven  to  not  provide  a  realistic  coupling  between  the  fluid  and 
the  structure,  other  codes  such  as  HRSR  (Hydraulic  Ram  Structural 
Response)  [28],  ERAM  or  EHRSR  were  developed  [29],  but  all  of 
them  showed  limitations  on  simulating  an  HRAM  event  since  none 
of  them  fully  coupled  the  mechanisms  of  fluid-structure  interac¬ 
tion  and  were  limited  to  simple  structures.  The  complicated  phys¬ 
ics  and  mechanics  of  HRAM  phenomena  were  not  satisfactory 
solved  until  higher-order  numerical  algorithms  were  incorporated 
into  the  codes  in  the  late  1980s.  Coupled  Euler-Lagrange  (CEL) 
methods  have  been  under  development  since  the  early-to-mid 
1990s.  They  combine  the  desirable  characteristics  of  Lagrangian 
and  Eulerian  formulations  so  that  the  distortion  problems  in  the 
fluid  are  avoided  due  to  the  Eulerian  approach  as  well  as  the  limits 
of  the  structure  are  perfectly  reproduced  by  means  of  the  Lagrang¬ 
ian  formulation.  These  methods  are  used  in  multiple  industries  for 
a  wide  variety  of  analysis  in  which  fluids  interact  with  structures 
or  when  high  distortions  may  appear  [21,30-37],  including  airbag 
and  tire-water  dynamics  in  the  automobile  field  [38,39],  the 


impact  of  bird  strikes  on  aircraft  [40,41  ],  and  the  effects  of  sloshing 
on  ships  [42]. 

In  the  last  years  and  motivated  by  different  industries  such  as 
aeronautics,  naval  or  more  recently  biomedical  sciences,  there 
have  been  new  advances  in  development  and  use  of  computational 
methods  for  fluid-structure  interactions  in  order  to  reach  more 
effective  computational  techniques  [43-48]  and  solving  more  dif¬ 
ficult  problems.  As  an  example  of  the  increasing  interest  on  solving 
industrial  fluid-structure  problems  the  works  of  Petitpas  et  al.  [49] 
or  Lecysyn  et  al.  [22,23]  can  be  mentioned,  in  which  a  ballistic  im¬ 
pact  on  an  industrial  tank,  filled  with  a  toxic  fluid,  is  studied. 

The  simulation  of  coupled  problems  of  fluid-structure  interac¬ 
tion  such  as  HRAM  added  to  the  modelling  of  carbon  fiber  compos¬ 
ites,  has  been  proven  to  be  a  complicated  task  and  is  still  quite 
challenging  [51,50].  The  suitability  and  predictive  capabilities  of 
different  techniques  (CEL,  ALE  or  SPH)  in  an  HRAM  problem  have 
not  yet  been  solved.  This  is  of  great  importance  since  analytical 
solutions  can  provide  only  a  limited  understanding  of  the  nature 
of  the  behavior. 

In  the  present  work,  numerical  simulations  of  a  water-filled 
CFRP  square  tube  subjected  to  impact  by  steel  spherical  projectiles 
at  different  velocities  are  shown.  The  simulations  are  performed 
with  the  commercial  code  ABAQUS/Explicit  V.6.10  using  a  Coupled 
Eulerian-Lagrangian  (CEL)  approach  to  reproduce  the  fluid-struc¬ 
ture  interaction  produced  as  a  consequence  of  the  HRAM  phenom¬ 
enon.  In  addition  a  woven  CFRP  material  model  is  implemented  by 
means  of  an  user  subroutine  to  appropriately  reproduce  the  behav¬ 
ior  of  the  structure.  Experimental  tests  data  regarding  pressure  in 
different  points  of  the  fluid,  cavity  evolution  for  different  impact 
velocities  and  the  failure  of  the  walls  are  compared  with  the 
numerical  results  in  order  to  assess  the  validity  and  accuracy  of 
the  proposed  numerical  model  in  reproducing  such  a  complex 
phenomenon. 

2.  Experimental  setup 

In  order  to  achieve  an  appropriate  and  wide  validation  of  the 
numerical  model,  the  same  authors  of  this  work  performed  exper¬ 
imental  tests  to  obtain  data  on  pressure,  cavity  evolution  and  fail¬ 
ure  of  the  CFRP  structure  due  to  the  HRAM  phenomenon.  The 
experimental  results  were  presented  in  a  previous  paper  [52]. 
These  tests  consisted  on  high  velocity  steel  sphere  impacts  against 
a  woven  CFRP  tube  filled  with  water.  In  Fig.  la  the  sketch  of  the 
experimental  device  used  in  the  mentioned  tests  can  be  seen. 

The  specimen  into  which  projectiles  were  impacted  corre¬ 
sponds  to  square  woven  CFRP  tubes  150  mm  wide,  2.2  mm  thick 
and  750  mm  long.  The  composite  woven  laminated  selected  was 
the  AGP-193-PW  manufactured  by  Hexcel  Composite,  composed 


(a) 


(b) 


Fig.  3.  (a)  Woven  CFRP  and  projectile  lagrangian  mesh  (b)  Impacted  zone  mesh. 
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Table  2 

Steel  and  PMMA  elastic  properties. 


P  (kg/m3) 

E  (GPa) 

V 

Steel 

7850 

210 

0.3 

PMMA 

1180 

3 

0.35 

by  10  plies  ([0]iO).  Each  ply  is  made  with  a  plain  weave  of  ASA 
fibers  and  the  8552  resin.  The  tube  was  closed  with  two  PMMA 
windows,  30  mm  thick,  fixed  with  four  steel  bars;  through  these 
windows,  the  whole  impact  and  penetration  process  was  recorded 
by  means  of  a  Photron  Ultima  APX-RS  digital  high-speed  camera;  a 
similar  setup  was  proposed  by  Nishida  et  al.  [53].  To  obtain  opti¬ 
mal  images  of  the  penetration  process,  it  is  necessary  an  appropri¬ 
ate  lighting,  which  was  provided  by  an  Arrisun  12  Plus  lamphead 
with  a  1200  W  Hydrargyrum  Medium  Arc  Iodide  (HMI)  lamp. 

In  order  to  obtain  pressure  values  inside  the  fluid,  two  pressure 
transducer  (PCB  138A06)  were  located  inside  the  tube.  The  posi¬ 
tion  of  the  pressure  gauges,  near  from  the  impact  point  (PTn) 
and  far  from  the  impact  point  (PTf),  can  be  seen  in  Fig.  lb.  In  adit- 
tion  six  uniaxial  strain  gauges  (350  Q,  2.120  Gage  factor  from 
Vishay  Measurements  Group  Inc.)  located  at  different  points  of 
the  CFRP  walls  can  also  be  seen  in  the  Fig.  lb. 

The  projectile  that  impacts  into  the  CFRP  tubes  is  a  12.5  mm 
diameter  steel  sphere;  it  is  accelerated  with  a  one  stage  light  gas 
gun,  capable  of  storing  gas  at  a  maximum  pressure  of  300  bar. 
The  length  of  the  barrel  is  4.5  m  and  its  calibre  25  mm.  The  speci¬ 
men  impacted  was  placed  inside  an  armored  steel  chamber  box 
lxlxl  m3.  The  chamber  had  a  small  circular  window  in  the 
front  for  the  projectile  to  pass  through,  and  two  large  lateral  win¬ 
dows  to  illuminate  the  specimen  and  capture  the  video  sequence 
of  the  impact.  Two  different  impact  velocities  were  performed: 
600  and  900  m/s.  Further  information  of  the  experimental  tests 
and  a  complete  detail  of  the  results  obtained  can  be  found  in  [52]. 

3.  Numerical  implementation 


problems  (drop  in  explicit  time  step,  worsening  in  results  accuracy, 
error  termination  of  simulations. . .  ).  By  the  other  hand,  an  Euleri- 
an  description  allows  the  material  to  flow  within  a  fixed  mesh, 
avoiding  mesh  distortion  problems.  Therefore,  in  this  work  the 
fluid  is  modeled  using  an  Eulerian  description  while  a  lagrangian 
description  is  used  to  model  the  surrounding  structure  and  the 
projectile. 

The  drawback  of  this  description  is  the  presence  of  convection 
terms  in  the  Eulerian  equations  [54].  The  operator  splitting 
proposed  by  Benson  is  a  very  efficient  method  for  solving  the 
Eulerian  equations  and  it  is  implemented  in  the  majority  of  the 
commercial  finite  element  codes  (also  in  ABAQUS).  In  this  meth¬ 
od,  the  step  is  divided  into  one  Lagrangian  step  and  one  Eulerian. 
As  a  consequence  of  this,  computational  costs  are  higher  than 
when  only  a  lagrangian  description  is  used.  The  coupling  algo¬ 
rithm  between  both  descriptions  is  performed  by  a  penalty  based 
contact  algorithm  [55].  Also,  it  must  be  taken  into  account  that  in 
order  to  let  the  fluid  move  through  the  Eulerian  domain,  the 
mesh  must  be  shared  between  the  fluid  and  other  material  or 
at  least,  void.  As  the  Eulerian  mesh  is  fixed  in  space,  if  only  one 
material  fills  that  mesh,  no  material  movement  will  be  seen  in 
the  simulations. 

3 A.  Carbon  fiber  woven  epoxy  laminates  model 

The  behavior  of  the  carbon/epoxy  woven  laminate  has  been 
modeled  as  an  orthotropic  elastic  material  until  failure.  This  kind 
of  approach  has  been  widely  used  in  impact  problems  on  com¬ 
posite  materials;  some  examples  are  the  Hou  et  al.  [56]  model 
for  tape  laminates  or  the  J.  Lopez-Puente  et  al.  [13]  model  for 
woven  laminates.  To  model  how  the  material  fails,  different 
damage  mechanism  are  defined:  fiber  failure,  matrix  failure 
and  delamination.  For  each  damage  mechanism  one  or  two  sca¬ 
lars  are  defined  by  means  of  the  actual  stress  tensor  and  the 
laminate  strength  properties;  its  value  could  vary  from  0,  which 
means  no  damage  to  1  which  means  fully  damage.  The  mecha¬ 
nisms  are: 


The  HRAM  phenomenon  is  a  highly  non-linear  transient  dy¬ 
namic  problem,  therefore  the  usage  of  an  explicit  finite  element 
code  is  recommended  to  try  to  reproduce  it.  In  this  work,  the  com¬ 
mercial  code  Abaqus/Explicit  v6.10  has  been  used  to  model  the 
aforemention  phenomenon  and  its  effects  on  a  CFRP  tube.  The  nat¬ 
ure  of  this  problem  makes  really  complex  the  employment  of  a 
lagrangian  method  in  which  the  finite  element  mesh  gets  de¬ 
formed  jointly  with  the  material.  In  situations  in  which  deforma¬ 
tion  is  extremely  large,  mesh  gets  distorted  leading  to  numerical 


The  fiber  failure  is  described  by  means  of  df i  and  dp,  one  for 
each  fiber  direction: 


dn  = 


if  (Tii  >  0 

^  if  ou  <  0 


(1) 


dp  = 


if  022  >  o 
if  022  <  o 


(2) 


97.2  mm 


Fig.  4.  Eulerian  domain  modeled. 
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Fig.  5.  Eulerian  mesh. 


Zc  •  @33  t  @33 

4  S13S23  Zc 


+  max 


(4) 


where  <7y  are  components  of  the  stress  tensor,  Si2,  Si3  and  S23  are 
the  shear  strengths  in  the  three  different  planes  and  finally  Zc  is 
the  strength  in  the  through-thickness  direction  under  compression. 
The  Eq.  (4)  applies  only  when  o33  <  0.  Eq.  (4)  is  a  modified  version  of 
the  expression  proposed  by  Hou  et  al.  [56]  for  this  mechanism  for 
tape  laminates. 

•  For  the  delamination  failure  mechanism  one  scalar  is  defined; 
the  expression  used  is  the  same  proposed  by  Hou  et  al. 


Table  3 

Water  properties. 


Property 

Unit 

Value 

Po 

(kg/m3) 

1000 

c 

(m/s) 

1448 

s 

(-) 

1.972 

V 

(Pas) 

0.00089 

r0 

(-) 

0.11 

P-m 

(J/kg) 

920.9 

where  <7n  and  cr22  are  the  stresses  in  the  warp  and  fill  direction 
respectively,  Xt  and  Xc  are  the  strengths  of  the  composite  laminate 
in  tension  and  compression  for  the  warp  direction,  and  finally  Yt 
and  Yc  are  the  strengths  in  tension  and  compression  for  the  fill 
direction. 

•  For  the  matrix  failure  damage  mechanism  two  parameters  are 
defined,  one  in  plane  direction  (dmi2),  an  the  other  one  in  the 
through-thickness  direction  (dm3).  The  corresponding  equations 
are: 


(3) 


where  again  Gy  are  components  of  the  stress  tensor,  and  Zr  is  the 
strength  in  the  through-thickness  direction  under  tension.  This  fail¬ 
ure  mechanism  is  computed  only  when  a33  >  0. 


Fig.  2  shows  the  stresses  that  participate  in  each  failure  mech¬ 
anism.  When  the  value  of  one  of  the  damage  parameters  described 
reaches  the  value  of  1,  the  components  of  the  stress  tensor  Oy  in¬ 
volved  in  the  failure  definition  are  set  to  zero.  To  avoid  numerical 
instabilities  during  the  simulations  due  to  sudden  changes  in  stiff¬ 
ness,  a  smooth  transition  is  used.  Eq.  (6)  shows  how  the  stress 
components  are  modified  as  a  function  of  the  corresponding  dam¬ 
age  parameters  defined.  For  the  numerical  simulations  accom¬ 
plished  in  this  work,  the  value  s  =  45  was  adopted. 


=  °ij 


2  -  es(du_1/2)\ 

2  -  es/2  ) 
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The  element  erosion  is  controlled  by  the  total  strain;  after  each  time 
increment  the  strain  tensor  is  calculated,  if  one  of  the  components 
reaches  a  critical  value,  then  the  element  is  removed.  The  properties 
of  the  carbon  fiber  woven  laminates  are  presented  in  Table  1. 


Fig.  6.  Impact  at  900  m/s.  (a)  Pressure  field  at  20  ps  after  the  impact,  (b)  Pressure  field  at  120  ps  after  the  impact,  (c)  Maximum  cavity  size  at  1.2  ms  after  the  impact,  (d) 
Contour  plot  of  stresses  in  fiber  direction  (<7n)  at  240  ps  after  the  impact. 
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Fig.  7.  Projectile  position  obtained  analytically,  numerically  and  experimentally  at 
900  and  600  m/s. 

3.2.  Box  and  projectile  lagrangian  FE  model 

The  model  of  the  problem  under  consideration  can  be  simplified 
attending  to  its  symmetry,  so  that  only  a  quarter  of  the  problem 
has  to  be  modeled  obtaining  a  desirable  reduction  of  the  computa¬ 
tional  cost  (Fig.  3a).  The  solid  parts  involved  in  the  problem,  the 
projectile  and  the  woven  CFRP  tube  which  is  closed  by  the  PMMA 
window,  are  modeled  as  follows. 

•  Woven  CFRP  tube.  It  is  discretised  by  means  of  eight  node  linear 
solid  elements  with  reduced  integration  and  hourglass  control, 
called  C3D8R.  Solid  elements  must  be  used  in  order  to  repro¬ 
duce  the  perpendicular  impact  and  penetration  process  into 
the  woven  CFRP  tube.  Mesh  density  is  higher  in  the  impacted 
zone  (lxl  mm2)  and  it  decreases  with  the  distance  from  the 
impacted  zone  (until  7.3  x  4.5  mm2),  performing  a  mesh  that 
accurately  simulates  the  damage  induced  according  to  previous 
works  [57].  A  detailed  of  this  refined  zone  can  be  seen  in  Fig.  3b. 
The  CFRP  tube  walls  present  five  elements  through  the  thick¬ 
ness,  so  each  element  models  two  plies.  The  material  model 
previously  described  in  Section  3.1  has  been  implemented 
through  a  user  subroutine  in  order  to  reproduce  the  woven 
CFRP  behavior. 


(a) 


•  Steel  projectile  and  PMMA  window.  The  projectile  is  discretised 
by  means  of  C3D8R,  while  the  PMMA  window  is  discretised 
by  means  of  four  node  conventional  shell  elements  with 
reduced  integration  (S4R).  The  element  size  of  the  projectile, 
although  relatively  bigger  than  those  elements  of  the  tube  in 
which  impacts,  allows  to  solve  in  a  properly  way  the  contacts 
with  the  tank  walls.  Both,  steel  and  PMMA,  are  modeled  as  elas¬ 
tic  materials  since  no  plastic  deformation  nor  damage  is 
observed  in  none  of  them  in  the  experimental  tests.  The  mate¬ 
rial  properties  used  in  the  model  for  PMMA  and  steel  can  be 
seen  in  Table  2. 

Finally  46,330  C3D8R  elements  have  been  used  for  the  woven 
CFRP  tube,  400  C3D8R  elements  for  the  steel  projectile  and  1968 
S4R  elements  for  the  PMMA  window. 

3.3.  Eulerian  fluid  model 

As  it  has  been  already  mentioned,  fluid  is  modeled  using  an 
Eulerian  description,  avoiding  numerical  instabilities  that  would 
have  appeared  using  a  lagrangian  approach. 

The  Eulerian  domain  used  in  the  simulations  can  be  seen  in 
Fig.  4.  The  domain  dimensions  have  been  defined  after  several 
studies,  so  that  the  deformed  CFRP  tube  walls  never  reach  the 
boundary  of  the  Eulerian  domain.  Therefore,  the  interaction  be¬ 
tween  fluid  and  the  surrounding  structure  is  computed  during 
the  whole  simulation  being  able  of  reproducing  the  effect  of  the 
HRAM  phenomenon. 

In  order  to  allow  the  fluid  inside  the  tube  to  flow  around  and 
interact  with  the  structure,  it  is  necessary  to  define  a  void  mesh 
outside  of  the  water  mesh  as  it  was  already  mentioned,  Fig.  4. 

The  whole  Eulerian  mesh  is  discretised  by  means  of  8  node  sol¬ 
ids  elements,  multimaterial,  with  reduced  integration  and  hour¬ 
glass  control,  called  EC3D8R. 

In  order  to  assure  an  accurate  contact  between  fluid  and  the 
surface  of  the  woven  CFRP  tube  walls,  both  parts  were  discretized 
with  the  same  element  size.  This  fact  avoids  possible  leakage  prob¬ 
lems  [58].  The  element  size  of  the  fluid  mesh  was  chosen  according 
to  previous  studies  regarding  deceleration  of  an  sphere  inside  a 
fluid.  A  simplified  model  without  the  CFRP  tube  was  used  to  per¬ 
form  several  impacts  and  compare  the  results  of  deceleration  of 
the  sphere  for  different  element  sizes  with  the  theoretical  case. 
The  value  of  2  mm  have  been  chosen  because  it  reproduces  accu¬ 
rately  the  deceleration  with  a  reasonable  computational  cost. 


Time  [ms] 


(b) 


Fig.  8.  Pressure  time  history  near  (PTn)  and  far  (PTf)  from  impact  point  in  a  tube  completely  filled  (a)  impacted  at  900  m/s  (b)  impacted  at  600  m/s. 


292 


J.A.  Artero-Guerrero  et  al.  /  Composite  Structures  96  (2013)  286-297 


Fig.  9.  Sequence  of  projectile  penetration  into  the  tube  at  600  m/s  and  900  m/s.  Images  taken  from  experiment  at  (a)  83  ps,  (b)  413  ps,  (c)  663  ps  and  from  numerical 
simulation  at  (d)  80  ps,  (e)  420  ps  and  (f)  660  ps.  The  images  from  the  experimental  impact  at  900  m/s  were  taken  at  (g)  83  ps,  (h)  413  ps  and  (i)  663  ps  and  the  numerical 
impact  at  (j)  83  ps,  (k)  413  ps  and  (1)  663  ps. 


Finally,  316,734  Eulerian  elements  have  been  used  in  the  model, 
Fig.  5. 

As  it  was  already  mentioned,  the  interaction  between  the  fluid 
and  the  structure  is  made  through  a  coupling  algorithm  based  on  a 
penalty  method.  In  Abaqus/Explicit  this  contact  is  used  when  the 
general  contact  option  is  active. 

When  one  material  penetrates  into  another,  the  penalty  meth¬ 
ods  acts  introducing  the  effect  of  a  spring  that  moves  the  interfaces 
to  a  non-penetrating  situation  [55]. 

The  fluid  material  behavior  is  defined  by  the  following  viscous 
constitutive  equation  [58]: 

G  =  2  rjs'  -PI  (7) 

where  y\  is  the  dynamic  viscosity,  s'  is  the  deviatoric  strain  rate,  P  is 
the  pressure  and  I,  the  identity  tensor.  The  pressure  P  is  related  with 
density  p  using  a  Mie-Griineisen  equation  of  state,  where: 

p=^(,-T2)+r«£-  «8» 

for  compressed  materials  and 

0) 


for  expanded  materials,  avoiding  negative  pressure  in  the  fluid.  In 
Eq.  (8),  rj  is  the  nominal  volumetric  compressive  strain  rj  =  1  -  p0/ 
p,  p  is  the  density  and  p0  is  the  initial  density;  c  is  the  speed  of 
sound  in  water,  and  s  is  the  slope  of  the  us  -  up  curve  where  us  is 
the  shock  velocity  and  up  is  the  particle  velocity.  r0  is  the  Griineisen 
gamma  and  Em  is  the  internal  energy  per  unit  of  mass.  Water  prop¬ 
erties  are  obtained  from  [59],  see  Table  3. 

4.  Results 

In  order  to  validate  the  model,  the  experimental  data  of  a  com¬ 
pletely  filled  tube  impacted  at  900  and  600  m/s,  obtained  by  Varas 
et  al.  [52],  will  be  compared  with  the  numerical  results.  The  capa¬ 
bility  of  the  CEL  approach  in  Abaqus/Explicit  to  reproduce  fluid- 
structure  interaction  problems,  as  well  as  the  HRAM  phenomenon, 
will  be  analyzed.  In  addition,  the  results  regarding  the  effects  in  the 
CFRP  structure  will  determine  if  the  proposed  material  model  is 
capable  of  reproducing  the  appropriate  behavior  of  the  tubes.  Six 
validation  variables  have  been  identified:  the  qualitative  reproduc¬ 
tion  of  the  HRAM  stages,  projectile  position  inside  the  tank,  fluid 
pressure  field,  cavity  evolution,  CFRP  strain  values  and  final  failure 
of  tanks. 


P  =  0 
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Fig.  10.  Detail  of  the  strain  time  history  of  G1  obtained  numerically  and 
experimentally  at  900  and  600  m/s. 


Time  [ms] 

Fig.  11.  Strain  time  history  of  G2  obtained  numerically  and  experimentally  at  900 
and  600  m/s. 

•  HRAM  stages.  Fig.  6  shows  how  the  numerical  model  qualita¬ 
tively  reproduces  the  four  well  known  HRAM  stages. 

-  Shock  phase.  The  event  that  characterizes  this  stage  is  the 
hemispherical  pressure  wave  produced  due  to  the  projectile 
impact  and  that  is  transmitted  through  the  fluid  to  the  whole 
structure.  The  mentioned  pressure  wave  at  20  ps  after  the 
impact  can  be  seen  in  Fig.  6a. 

-  Drag  phase.  The  projectile,  while  travels  through  the  tank, 
transmits  part  of  its  kinetic  energy  to  the  fluid;  hence  the 
fluid  is  displaced  from  the  projectile  path  generating  a  radial 
pressure  field  and  a  cavity  in  the  projectile  wake,  Fig.  6b. 

-  Cavity  phase.  This  stage  is  represented  by  the  growth  of  the 
cavity  generated  by  the  projectile.  Fig.  6c  depicts  the  maxi¬ 
mum  size  reached  by  the  cavity.  The  cavity  pushes  the  fluid 


Fig.  12.  Strain  time  history  of  G5  obtained  numerically  and  experimentally  at  900 
and  600  m/s. 


against  tube  walls,  inducing  significant  deformations  in  the 
structure.  This  stage  is  the  main  responsible  for  the  most 
important  damage  effects  on  the  structure.  This  is  very  well 
reproduced  by  the  numerical  model. 

-  Exit  phase.  The  exit  wall  of  the  tank  is  pre-stressed  before  the 
projectile  go  through  it.  This  is  due  to  the  initial  shock  stage 
and  the  subsequent  loading  by  the  fluid.  Fig.  6d  shows  how 
the  model  reproduces  this  behavior.  The  exit  wall  is  already 
stressed  before  the  projectile  impact,  which  would  explain 
why  damage  and  strains  are  larger  in  the  exit  wall  than  in 
the  entry  one. 

•  Projectile  position.  As  it  was  already  mentioned,  the  experimen¬ 
tal  projectile  position  can  be  obtained  by  means  of  the  high¬ 
speed  camera.  This  data  have  been  compared  with  numerical 
and  analytical  projectile  position,  obtained  integrating  from 
Newton’s  second  law  [60] 

mP^  =  -lPwAoCdV2p  (10) 

where  mp  and  Vp  are  the  projectile  mass  and  velocity,  pw  is  the  fluid 
density,  A0  is  the  projected  frontal  area  of  the  projectile  and  Cd  is  a 
dimensionless  drag  coefficient.  According  to  the  range  of  velocities 
considered,  a  value  of  0.4  for  Cd  was  chosen. 

The  experimental  and  numerical  time  history  of  the  projectile  is  de¬ 
picted  in  Fig.  7.  A  good  correlation  between  the  curves  it  is  ob¬ 
served,  showing  a  maximum  difference  of  about  a  6%.  The  trend 
of  the  projectile  velocity  can  be  observed  by  means  of  the  projectile 
position  slope  in  Fig.  7.  The  velocity  decreases  inside  the  fluid, 
transforming  part  of  its  kinetic  energy  into  pressure  and  kinetic 
energies  in  the  fluid.  As  it  was  mentioned  before,  the  projectile  po¬ 
sition  correlation  was  used  to  determine  the  mesh  discretization  in 
the  projectile  direction.  It  has  to  be  noted  that  with  a  finer  mesh,  as 
other  authors  observed  [58,61,62],  better  results  can  be  obtained. 
Nevertheless  instabilities  in  the  contacts  between  the  fluid  and 
the  structure  as  well  as  the  computational  cost  growth,  would  made 
useless  the  improvement. 

•  Pressure  field.  As  was  already  mentioned,  two  pressure  trans¬ 
ducers  were  located  inside  the  fluid  at  different  positions;  this 
allows  knowing  how  the  distance  from  the  impact  point  affects 
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Fig.  13.  Final  failure  at  the  exit  wall  of  the  CFRP  tube.  v=  900  m/s. 
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Fig.  14.  Final  failure  at  the  exit  wall  of  the  CFRP  tube.  v=  600  m/s. 
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Fig.  15.  Detail  of  the  final  failure  of  the  CFRP  tube  at  900  m/s. 


(a)  (b) 


Fig.  16.  a)  Maximum  pressure  in  impacts  at  450  m/s,  600  m/s,  750  m/s  and  900  m/s.  (b)  Maximum  cavitation  size  in  impacts  at  450  m/s,  600  m/s,  750  m/s  and  900  m/s. 
Obtained  from  numerical  results. 


the  pressure  generated  in  the  fluid.  One  pressure  transducer 
(PTn)  was  located  at  30  mm  from  the  wall  and  75  mm  away 
from  the  shot  line,  that  is  81  mm  from  the  impact  point.  The 
other  transducer  (PTf)  was  placed  in  the  middle  of  the  tube, 
150  mm  away  from  the  projectile  trajectory  and  171  mm  from 
the  impact  point,  see  Fig.  lb. 

The  pressure  values  obtained  in  the  pressure  transducers,  both 
experimental  and  numerically,  when  a  completely  filled  tube  is 
impacted  at  different  velocities  are  depicted  in  Fig.  8.  A  good 
correlation  can  be  observed  not  only  in  peak  values  but  in  the 
trends  regarding  impact  velocity  and  pressure  transducer 


location  influences.  It  can  be  observed  how  the  intensity  of 
the  pressure  wave  increases  with  the  impact  velocity,  being 
the  pressure  peak  value  at  900  m/s  double  than  at  600  m/s;  this 
result  is  related  to  the  fact  that  the  pressure  is  proportional  to 
kinetic  energy  [24].  Pressure  values  in  the  far  pressure  trans¬ 
ducer  are  smaller  because  the  energy  of  an  spherical  pressure 
wave  diminishes  with  the  distance  to  its  source.  Nevertheless 
it  has  to  be  noted  that  better  predictions  are  achieved  in  points 
near  impact  (PTn)  than  far  from  the  impact  (PTf);  this  is  proba¬ 
bly  due  to  the  fact  that  the  mesh  is  finer  near  the  impact  point. 
As  pressures  are  averaged  over  the  elements,  the  usage  of  larger 
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Velocity  [m/s] 

Fig.  17.  Size  of  horizontal  and  vertical  cracks  appeared  in  a  quarter  of  the  exit  wall 
in  impacts  at  450  m/s,  600  m/s,  750  m/s  and  900  m/s.  Obtained  from  numerical 
results. 

elements  leads  to  lower  values  of  pressure.  Finally,  the  begin¬ 
ning  of  the  pressure  pulse  is  also  well  predicted  by  the  numer¬ 
ical  model. 

•  Cavity  evolution.  Since  cavity  is  the  main  cause  of  the  deforma¬ 
tion  and  failure  of  the  tanks,  it  is  interesting  to  analyze  how  the 
cavity  evolution  is  predicted  by  the  numerical  model.  Numeri¬ 
cal  cavity  evolution  can  be  compared  with  the  images  of  the 
penetration  process  recorded  with  the  High-Speed  Camera.  In 
Fig.  9,  experimental  and  numerical  images  at  different  impact 
velocities  are  shown.  The  images  show  the  initiation  and  the 
subsequent  growth  of  the  cavity  inside  the  fluid.  Good  correla¬ 
tion  between  experimental  and  numerical  cavity  evolution  is 
observed  for  both  impact  velocities,  Fig.  9.  Also,  it  can  be  seen 
that  the  cavity  grows  faster  when  the  impact  is  at  900  m/s, 
therefore  the  energy  of  the  fluid  is  higher  and  the  damage 
effects  in  the  structure  are  more  severe.  According  to  that, 
deformation  and  final  failure  of  the  CFRP  tank  are  influenced 
by  the  projectile  velocity. 

•  CFRP  strain  values.  As  it  has  been  said  before,  six  uniaxial  gauges 
were  attached  to  the  exterior  surface  of  the  CFRP  tube.  Two 
were  located  in  the  entry  wall:  G1  and  G2  located  at  50  mm 
and  1 50  mm  from  the  impact  point  respectively.  Other  two 
strain  gauges,  G3  and  G4,  were  located  in  the  exit  wall,  in  the 
same  positions  as  G1  and  G2.  Finally  two  more  gauges,  G5 
and  G6,  were  located  in  the  middle  of  the  lower  and  upper  wall, 
respectively  (Fig.  lb).  A  comparison  between  the  numerical  and 
experimental  strain  gauges  results  are  shown  in  Figs.  10-12.  It 
has  to  be  mentioned  that  some  of  the  experimental  signals 
are  not  longer  valid  few  time  after  the  impact.  This  is  due  to 


the  violent  shake  of  the  tube  when  it  is  impacted  and/or  the 
peeling  of  the  outer  plies,  which  may  induced  the  debonding 
of  the  gauges.  In  addition,  it  is  worth  to  say  that  the  positive 
sign  in  strain  data  indicates  that  the  wall  is  displaced  outwards 
Fig.  10  depicts  the  strain  data  at  gauge  G1  (closest  gauge  to  the 
impact  point)  for  a  completely  filled  tube  impacted  at  900  and 
600  m/s.  Although  the  experimental  data  is  no  longer  valid 
few  time  after  impact,  it  can  be  observed  that  the  numerical 
model  reproduces  the  trend  observed  experimentally  in  the  first 
instants  of  time.  The  strain  data  of  gauge  G2,  placed  in  the  entry 
wall  far  from  the  impact  point,  is  shown  in  Fig.  1 1.  It  can  be  said 
that  the  numerical  results  are  in  accordance  with  the  experi¬ 
mental  data,  both  in  maximum  values  and  trends.  Finally  the 
strain  time  history  in  the  middle  of  the  upper  wall  and  directly 
above  the  trajectory  of  the  projectile  (G5)  is  shown  in  Fig.  12.  As 
in  the  other  gauges,  the  numerical  predictions  agree  with  the 
trends  observed  experimentally.  It  has  to  be  noted  that  the 
experimental  curve  obtained  at  900  m/s,  diminishes  few  time 
after  impact.  This  is  probably  due  the  effect  of  the  crack  that 
appears  in  the  edge  which  provokes  the  relaxation  of  the  wall 
as  is  pointed  in  [52].  Although  the  crack  also  appears  in  the 
numerical  simulation,  it  does  not  affect  the  whole  thickness  of 
the  tube  wall,  hence  it  is  not  deep  enough  to  isolate  the  wall 
from  the  behavior  of  the  rest  of  the  tube.  Finally  it  is  worth  to 
mention  that  the  instant  in  which  the  maximum  strain  values 
appears  in  Figs.  11  and  12  matches  with  the  time  in  which 
the  maximum  cavity  size  is  reached  (1.2  ms  after  the  impact). 
This  confirms  the  importance  of  the  effect  of  the  cavity  in  the 
behavior  of  the  structure,  being  more  important  than  the  pres¬ 
sure  pulse. 

•  Failure  in  tubes.  Due  to  the  complexity  of  this  process  and  the 
influence  of  CFRP  manufacturing  in  the  tube  behavior,  it  is  not 
easy  to  reproduce  same  experimental  failure  in  numerical  sim¬ 
ulation.  Nevertheless,  similar  trends  can  be  seen  in  both  exper¬ 
imental  and  numerical  final  failures.  It  can  be  seen  in  Fig.  13 
that  the  most  damaged  area  is  the  exit  wall,  where  a  cross-type 
failure  appears.  When  the  tube  is  impacted  at  a  lower  velocity, 
600  m/s,  the  mentioned  cross-type  failure  also  appears  in  the 
exit  wall  with  a  smaller  size,  Fig.  14.  This  behavior  can  be  seen 
both  in  the  experimental  and  numerical  impacts,  although  the 
numerical  length  is  underestimated.  Fig.  15  shows  how  the 
numerical  model  predicts  the  longitudinal  crack  that  appears 
on  the  edge  of  the  tube  in  experiments.  In  this  case  the  model 
overestimates  the  length  of  such  a  crack. 

5.  Influence  of  velocity  in  HRAM  effect 

Once  the  numerical  model  has  been  validated  with  experimen¬ 
tal  results,  the  influence  of  velocity  in  HRAM  phenomenon  is  ana¬ 
lyzed  using  numerical  simulations.  Therefore,  four  numerical 


(a)  (b) 

(C)  (d) 


Fig.  18.  Horizontal  and  vertical  cracks,  represented  in  blue,  appeared  in  a  quarter  of  the  exit  wall,  (a)  450  m/s,  (b)  600  m/s,  (c)  750  m/s,  (d)  900  m/s.  Obtained  from  numerical 
results.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 
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impacts  have  been  simulated  at  different  projectile  velocities  (450, 
600,  750  and  900  m/s).  Maximum  peak  pressure  wave,  maximum 
size  of  cavity  and  failure  in  tubes  are  compared  between  the  differ¬ 
ent  impacts.  Fig.  16a  shows  the  influence  of  the  impact  velocity  on 
the  maximum  peak  pressures,  so  that  when  the  velocity  is  doubled, 
the  pressure  is  approximately  four  times  higher.  Therefore,  the 
pressure  wave  magnitude  is  proportional  to  the  projectile  kinetic 
energy  as  Lecysin  et  al.  mentioned  in  [24].  In  far  pressure  trans¬ 
ducer  the  peak  pressure  is  smaller,  as  it  was  already  mentioned, 
because  of  the  higher  distance  to  the  wave  source.  Fig.  16b  shows 
the  relation  between  impact  velocity  and  maximum  cavity  size.  A 
linear  relation  between  them  can  be  observed  until  an  impact 
velocity  of  750  m/s  is  reached.  Above  this  velocity,  the  slope  is  re¬ 
duced  due  to  the  fact  that  the  cavity  is  about  the  size  of  the  tube. 
Therefore,  when  the  impact  velocity  increases,  the  pressure  and 
the  cavity  size  increase  causing  a  higher  structural  damage,  as  it 
can  be  seen  in  Figs.  17  and  18.  Fig.  17  depicts  the  relation  between 
the  extent  of  vertical  and  horizontal  cracks,  caused  by  fiber  failure 
in  both  in-plane  axes,  and  the  impact  velocity.  In  Fig.  18,  it  can  be 
seen  the  failures  appeared  in  a  quarter  of  the  CFRP  tube.  No  differ¬ 
ent  failure  mechanism  can  be  observed  in  the  four  images, 
although  a  failure  in  the  curved  edge  is  observed  at  900  m/s, 
Fig.  15. 

6.  Conclusions 

The  capability  of  the  used  numerical  technique,  CEL  approach, 
in  modelling  the  fluid  structure  interaction  was  analyzed,  as  well 
as  the  accuracy  of  the  CFRP  material  model  proposed,  in  reproduc¬ 
ing  the  experimental  failures  appeared  in  a  HRAM  event.  Numeri¬ 
cal  results  were  compared  with  the  experimental  ones  obtained  in 
a  previous  work  done  by  Varas  et  al.  [52].  In  addition,  various 
numerical  impacts  have  been  performed  to  analyze  the  influence 
of  impact  velocity  in  HRAM  effects.  Some  conclusions  have  to  be 
remarked: 

•  Coupled  Eulerian  Lagrangian  approach  implemented  in  ABA- 
QUS/Explicit  has  been  shown  as  a  reliable  tool  to  reproduce 
HRAM  phenomenon.  The  model  is  able  to  reproduce  HRAM 
stages:  shock,  drag,  cavity  and  exit  phases  from  qualitative 
and  quantitative  perspective.  The  projectile  position  is  well  cap¬ 
tured,  reproducing  the  energy  exchange  between  the  projectile 
and  the  fluid.  The  cavity  evolution,  main  cause  of  the  tank  final 
failure,  is  accurately  reproduced  as  compared  with  the  images 
taken  from  the  experimental  tests. 

•  Water  material  model  used  reproduces  accurately  the  peak  val¬ 
ues  generated  by  the  pressure  wave  in  the  fluid.  Also  it  has  been 
shown  how  the  magnitude  of  the  pressure  wave  diminishes 
with  the  distance  from  the  impact  point. 

•  The  CFRP  model,  implemented  by  an  user  subroutine,  repro¬ 
duces  the  experimental  failures  in  the  most  damaged  area  in 
the  tube.  It  has  been  shown  how  the  numerical  strain  values 
registered  agrees  with  the  trends  observed  experimentally. 
Also,  strain  measures  registered,  both  numerical  and  experi¬ 
mental,  shows  that  the  maximum  strain  values  occurs  when 
the  cavity  has  reached  its  maximum  size,  confirming  the  impor¬ 
tance  of  the  cavity  effect  in  the  behavior  of  the  structure  and 
hence  in  the  final  tank  failure. 

•  A  higher  velocity  impact  produce  higher  pressure  peaks  in 
shock  phase  and  also  a  larger  cavity.  Therefore,  more  damaging 
HRAM  effects,  as  it  is  shown  in  experimental  and  numerical 
tests.  The  relation  between  magnitude  of  pressure  pulse  and 
impact  velocity  is  quadratic,  while  maximum  cavity  size  is 
linear. 
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